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The Bose-Hubbard model in an external magnetic field is investigated with strong-coupling 
perturbation theory. The lowest-order secular equation leads to the problem of a charged particle 
moving on a lattice in the presence of a magnetic field, which was first treated by Hofstadter. We 
present phase diagrams for the two-dimensional square and triangular lattices, showing a change 
in shape of the phase lobes away from the well-known power-law behavior in zero magnetic field. 
Some qualitative agreement with experimental work on Josephson-junction arrays is found for the 
insulating phase behavior at small fields. 
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I. INTRODUCTION 

The simplest model of strongly interacting bosons is the Bose-Hubbard modeL (BHM) , which has been used to 
describe superfluid heliumtl, Cooper pairs in thin granular superconducting filmsa and Josephson iimction arraysa. 
Much theoretical work has concentrated on the phase diagram of the BHM in zero magnetic fieldel because of the 
technical problems associated with introducing an external magnetic field (such as the sign problem in quantum 
Monte Carlo simulations). Experimentalists on the other hand have concentrated on studying systems in an external 
magnetic field because the phase transition can be tuned by adjusting the magnetic field rather than changing the 
samples measurecoj. We study the BHM on two-dimensionaLlattices in a perpendicular magnetic field by extending 
the strong-coupling perturbation theory for the field-free cased. This theoretical technique can incorporate magnetic- 
field dependence in a straightforward manner and is useful in studying field-tuned transitions. We concentrate on 
pure systems in this contribution and do not include any effects due to disorder. We present zero-temperature phase 
diagrams, study excitation-gap energies, calculate the dynamical critical exponent zv for small fields, and compare 
our theoretical results with experimental ones. 

The BHM contains the key physics of a many-particle bosonic system with competing potential and kinetic energy 
effects. The typical zero-temperature phase diagram for the noijwmagnetic case shows incompressible Mott-insulating 
(MI) phases surrounded by compressible superfluid phases (SF)cl. The insulator to superfluid transitions at the tip of 
the lobes (where the density remains constant) are driven by quantum phase fluctuations, while those at the sides of 
the lobes (where the density varies) are driven by density fluctuations, i.e. particle or hole excitations. Introducing a 
magnetic field is expected to increase the region of the MI phase because the localizing effect on the itinerant bosons 
reduces the stability of the SF phase. 

We consider bosons with total spin 0. The only effect of a perpendicular magnetic field H is then on the orbital 
motion of the bosons, which effects changes in the phase of the hopping-matrix T = (tjk) between lattice sites j 
and fc. By choosing a Landau-gauge for the vector-potential A(r) = H(0,x,0), the Bose-Hubbard Hamiltonian in an 
external magnetic field becomes 

Hbhm = - ( f jk b ] b k + h - c -) + ^oinj - I) - V^rij, (1) 

<jk> j j 

where the hopping matrix is nonzero only between nearest neighbors and is given by 

t jk = te-^ A >» , A jk = -L / A(r) ■ dr. (2) 

00 Jj 

This hopping matrix is Hermitian because t is real and Ajk = —Akj. Here the boson creation operator for lattice site 
j is bp hj — bjbj is the corresponding number operator, U is the on-site repulsion of the bosons, and fi is the chemical 
potential. We choose U to be our energy scale and measure all energies in units of U. The magnetic flux quantum is 
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given by (f>o — — , and the magnetic flux per plaquette 2ncf) — § A(r) ■ dr ~ a 2 H is a measure of the strength of the 

magnetic field H (where a is the lattice spacing). _ 

The form of the zero-temperature phase diagram can be understood by starting from the atomic limitQ, where the 
hopping t is zero and every site is occupied by a fixed number of bosons no. The energy to add one boson onto a 
site with n bosons is E(n + 1) — E(n ) = n U — fi, so that there is a finite energy gap when (n — 1)U < fi < n U 
and the system is an incompressible Mott insulator. When fi = n U, then all states with a density between no and 
n + 1 bosons per site are degenerate in energy, and the system becomes a compressible fluid. As the strength of 
the hopping matrix elements increases, the range of the chemical potential about which the system is incompressible 
decreases. The Mott insulator phase disappears at a critical value of the hopping matrix elements (which depends on 
the strength of the magnetic field and the density of the insulating phase) and the system becomes a superfluid (see 
Fig.^j). We provide a systematic study of the BHM in a magnetic field by examining the system in a perturbative 
expansion about the atomic limit with the boson kinetic energy acting as the perturbation. 

The paper is organized as follows: Section II describes the formalism used in the strong-coupling perturbation 
theory in the presence of a magnetic field. Section III presents the results for the phase diagrams and the excitation 
energies in a magnetic field, and Section IV contains the conclusions. 



II. FORMALISM 



Our procedure is to calculate the ground-state energy of the MI phase with n bosons per site E g (jiQ,t), and of 
the excited states in the charge sector with one extra boson E p (no,t) and one extra hole Eh{riQ,t) 1 in a Rayleigh- 
Schrodinger perturbative expansion in the hopping matrix element t (the kinetic energy is chosen as the perturbative 
part of the Hamiltonian) . When the energy of the MI and the state with one extra boson are equal, the system 
undergoes a phase transition from the incompressible MI phase to the compressible SF phase with density larger than 
no. The similar occurs when the state with one extra hole is degenerate with the MI phase (except now the density of 
the SF phase is less than n«). The detailed formalism of the strong-coupling expansion for the ground-state energies 
has already been presented!! The only modifications of the previous calculations needed here are to take into account 
the fact that the hopping matrix now has a complex phase and the changes required for the nonbipartite hopping 
matrix of the triangular lattice. The important parameter that enters the results is the minimal eigenvalue et of the 
kinetic energy matrix — tjk which includes a factor e that depends on the magnetic field. This parameter determines 
how the degeneracy is lifted in the first-order secular equation for the energy of the excited states in the charge sectox. 
Formally this solution of the minimal eigenvalue is identical to finding the band minimum in the Hofstadter problemH. 

The Mott phase diagram is determined by the two Mott phase boundaries — one for the particle excitations [where 
E p (no, t) — E g (no, t) = 0] and one for hole excitations [where E^ijiQ, t) — E g (no, t) — 0]. For each value of t there is a 
critical value of the chemical potential where the system changes phase from an incompressible to compressible fluid. 
The upper and lower curves for the Mott phase lobe are then described by this critical value of the chemical potential 
H p /h{t)- The results of our expansion through third-order are 

fjtp(t) =n + e(n + l)t - n (n + l)e 2 t 2 + y (5n + A)zt 2 

25 7 

- n (n + l)[(2n + l)(-2 - e 2 ) + z(—n Q + -)]et 3 

+ 12<J tet n (y n 2 + yn + 3) cos(27T0)i 3 + (9(i 4 ) 

Hh{t) = n o - 1 - en> t + n (no + l)e t — (5n + l)zt 

25 11 

+ n (n + l)[(2n + l)(-2 - e 2 ) + z(—n Q + —)]et 3 

- 12* Iat (no + l)(^ng + 5n + i) cos(2^)t 3 + (9(t 4 ), (3) 

where Si a t = or 1 for the square or the triangular lattice, respectively, and z is the corresponding number of nearest 
neighbors (4 or 6, respectively). These results have been verified by both analytical small-cluster calculations and by 
numerical cluster expansions. 

The magnetic field appears in two places: (1) The magnetic field couples to the orbital current of the bosons as 
the particle or hole travels around a lattice plaquette and encloses the flux 2n(f>. In our third-order calculation, the 
orbital coupling only enters for the triangular lattice (it enters at fourth order for the square lattice as four hops are 
required to enclose a plaquette). (2) The other effect of the magnetic field is to change the minimal energy of the 
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FIG. 1. Band minimum of the magnetic band structure e(cj>) for the square and triangular lattice. 



extra particle or hole moving in the Mott phase background. Although the location of this mininum in the Brillouin 
zone is gauge-dependent, its value is not (thus e = mm^ eB e((fc))). 

We consider a rationaLilux 27rr/> = 2Ttp/q. To find the minimal eigenvalues et of the hopping matrix T, we follow 
Bellisardlj and HasegawaE2l In the Landau gauge, the system maintains its translational invariance in the x-direction 
while it requires m steps for translational invariance in the y-direction. For the square lattice, the period m equals q, 
for the triangular lattice, m — q/2 for even q and m — q for odd q. A Fourier-transformation now changes T to the 
following m x m- matrix T^(fc): 



f+{k) = t 



( Mi A 1 
A* M 2 ••• 

'■. '•. 



with 



M„ q = -2 cos(k y a + 27mr/»), 
M*" = -2cos(* y o + 47rn^), 











(4) 



A s i — 

A n — 

At" = 



Ak x a 



ik 



^i27r(f)(2n-\-l)-\-ik y a 



(5) 
(6) 



The eigenvalues within the corresponding Brillouin-zonc IB = {k | < k x < 2tt/iti, < k y < 2ir/m} determine the m 
energy bands of a boson moving in the magnetic field. The minimal eigenvalue of the band-structure et is shown in 
Fig. |l| as a function of the magnetic flux per plaquette 2n(j>. This is the parameter that enters Eq. (Q) to determine the 
Mott phase boundary. Notice how on the square lattice the largest dip is at = 1/2 followed by 1/3, 1/4, and 2/5, 
... while on the triangular lattice the sequence corresponds to 1/2, 1/4, 1/3, and 3/8, ... The relation e(r/>) = e(l — r/>) 
holds for all lattices, since the flux of 27r is equivalent to a flux of 0. Hence, the maximal magnetic field that can 
be applied corresponds to = 1/2. This maximal field configuration is realized by all real hopping matrix elements 
for each lattice: on the square lattice one takes three positive and one negative matrix element on each plaquette, 
while on the triangular lattice one takes all matrix elements to be negative. In particular, the "fully frustrated" case 
of <fi = l/2_pn a triangular lattice is the only nonzero magnetic-field case that can be easily treated by a high-order 
expansionEj because it maintains the full periodicity of the triangular lattice in zero magnetic field. 

Rather than plotting the phase diagrams for just a third-order expansion, we choose to extrapolate our results 
using knowledge about the overall structure of the phase diagram. There are numerous ways in which one can 
envision extrapolating the results of our third-order expansion to higher order. It has been demonstrated by Elstoer 
and Monien that a Pade analysis of the strong coupling perturbation series yields rapid convergence in ze*Q fieldE3. 
We apply this method to the magnetic-field case for small magnetic fields by using a Pade approximantEj for the 
logarithmic derivative of the difference in particle and hole Mott phase boundaries A(t) = (J, p (t) — ^h{t). We assume 
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FIG. 2. Phase diagram for the square (a) and the triangular (b) lattice in relatively small magnetic fields using the Pade 
analysis. 



the same behavior as found in the zero-field case, where the tip of the Mott lobe has a power law "critical point" with 
an exponent zv, A(t) = A(t)(t cr u — t) zv . Then the Pade analysis fits 

3, \ / \ zv A'it) 
at t - t crit A(t) 

to estimate the critical point and the dynamical critical exponent (the pole determines t cr n and the residue determines 
zv). The rest of the Pade approximant determines A(t), which then allows A(t) to be constructed. A similar Pade 
analysis for the midline of the Mott lobe fi m (t) = ^[Mp(^) + (which is a regular function of t) then allows the 

entire phase diagram to be constructed. 



III. RESULTS 

Atomic systems with a lattice spacing a of around 2A require a field of H = 5 x 10 3 Tesla for a half flux quantum 
per plaquette. This is too large a field to be accessed experimentally, hence atomic systems always lie in the low-field 
region, where the perturbation theory is most accurate. However for macroscopic lattice systems, with a = 2 x 10~ 4 cm 
(as in a Josephson-junction array) the whole magnetic field range is attainable by experiments with fields as low as 
H = 0.5 Gauss. We show below how .our theoretical results compare to the superfluid-insulator transitions on two 
dimensional Josephson-junction arraysO. 

Fig. H presents the phase diagram of the first three lobes for the square lattice §(a) and the triangular lattice ||(b) 
in the low magnetic-field region. The incompressible, Mott-insulating phase grows in size when the magnetic field 
increases from zero. This shows, as expected, the localizing effect of the magnetic field on the itinerant bosons. It 
appears that the shape of the lobes changes from a simple power-law dependence at the tip to something else. There 
are three possibilities for the shape of the phase lobe at the tip: (i) it remains power-law-like with an exponent 
that may depend on the magnetic field; (ii) it has a more exotic "cusp-like" dependence as seen in the Kosterlitz- 
Thouless transition; or (iii) it has a discontinuous change in slope at the tip, and thereby has a "first-order-like" 
shape corresponding to the crossing of the two curves representing the upper and lower phase lobes, with no "critical" 
behavior at the tip (this last result corresponds to zv = I). There is theoretical evidence in favor of this last conclusion, 
which is a [Similar behavior to what happens to the MI phase in the presence of disorder, but here the explanation 
is differenttil. If <fi = p/q is expressed in lowest terms, then the order parameter requires q components to describe 
it. As q increases, it is more likely that the transition is first-order rather than second-order, and hence one would 
expect that the tip has a slope discontinuity immediately upon the introduction of the magnetic field, since small 
fields correspond to large values of q. The Pade analysis given above assumes that (i) holds, but as zv approaches 
I for larger fields, the system crosses over to the behavior of (iii) unless something more exotic like (ii) intervenes. 
We cannot distinguish between scenarios (i) or (iii) with the low-order expansions presented here, but our results 
suggest scenario (iii) is correct, because the exponent zv rapidly approaches 1 as the field is increased. Clearly further 
theoretical analysis is needed to decide this issue. 

For larger magnetic fields, we find that the two boundaries of the Mott phase lobe no longer cross at a critical value 
of t, but rather they "repel" each other. This indicates that the tip of the MI phase lobe has moved out to such a large 
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FIG. 3. (a) Dependence of the gap energy (for fixed hopping) on the magnetic field. Notice the kinks that develop at 
commensurate magnetic fields (4> = 1/3 or 1/4). The curves are cut off at <f> m 0.35 because the perturbation theory fails for 
larger values of <j>. (b) Dependence of the dynamical critcal exponent zv on magnetic field when the truncated expansions, or 
a Pade approximant, are fit to the power-law behavior at the tip of the lobe. Notice that zv appears to be field-dependent and 
rapidly approaches 1. 



value of t that the perturbation theory needs to be carried out to higher order to achieve proper convergence. Such 
an exercise has been-already carried out through a linked-cluster expansion for the fully frustrated case (0 = 1/2) on 
a triangular latticellj. The calculation from 3rd to 11th order converges toward a "first-order" transition (zv = 1) at 
t cr n ~ 0.06 which yields a much larger insulating phase regime than with no magnetic field where t cr n ~ 0.038. A 
comparable experimental result is foundcl where the insulator to superfluid transitions on Josephson-junction arrays 
were studied by changing the ratio of the Josephson coupling energy Ej and the charging energy Eq. Since Ej is the 
energy scale for the superconducting coupling between the islands and Eq is the scale for the interaction between the 
charge carriers, they can be related to the hopping energy t and the on-site repulsion U, respectively. The results for 
small temperatures can be extrapolated to T = and show a transition at larger Ej/Ec for the 0=1/2 case than 
for the — case, which qualitatively matches with our calculations (since t cr it/U increases with magnetic field). 

The experimental data can be analyzed in one of two dual pictures with the Bose Hubbard model. In the first view, 
the bosons are the Cooper pairs on each island, and the main effect of the magnetic field is to modify the hopping 
matrix. Vortices appear in this case as supercurrent loops in the system. While the density of vortices increases as the 
magnetic field is modified, the density of Cooper pairs remains essentially constant over the low field ranges (on the 
order of Gauss) explored in the experiment^ The second picture employs the duality transformation that views the 
vortices themselves as the bosonic particlesuj. The mobility of the vortices is determined by Eq and their interaction 
by Ej , so the roles of those parameters are reversed in this case. The magnetic field then takes on the role of the 
chemical potential. We will not employ this second picture here, which has been used to evaluate the Mott insulating 
phase of the vortices in the quasi-one-dimensional Josephson junction arraystZI. 

In the following we compare small- field results. As shown in Fig. ||the MI phase area and the location of the tip oi 
the lobe increase with increasing magnetic field up to about = 0.1 (for the square lattice). The measurements ofEl 
show (Ej/ Ec)crit ~ 0.59 and 0.83 for = and 0.1 respectively, which is similar to our results with t cr u increasing 
by about 45% from 0.057 to 0.083 as increases from to 0.1. In addition, the tip location for larger field strengths 
saturates in the experiments, which is a result that we cannot confirm, because our analysis through third-order in t 
fails when t cr u becomes too large. 

In Fig. |^(a) we show the evolution of the excitation-gap energies A(t) for fixed hopping t as a function of magnetic 
field. Initially the gap energy increases which indicates the increasing stability of the MI phase in a weak magnetic 
field. For larger magnetic fields, we see dips in the gap energy around rational magnetic fluxes such as = 1/3 or 1/4 
(which can be reliably calculated for small values oft). A similar commensurable structure js found in experiments 
that measure the zero-bias resistance (Rq) of a Josephson-junction array in a magnetic fielcH. A small i?o indicates 
proximity to the SF phase. At rational the vortices form a lattice which is commensurable with the Josephson- 
junction array and favors a pinning of the vortices, leading to a decrease of the zero-bias resistance. Qualitatively, we 
see the same behavior in the excitation gap energies and the zero-bias resistance. Notice how the dips track closely 
with the dips seen in the minimal eigenvalue of the hopping matrix as shown in Fig. [l| This behavior explains the 
dips seen in the experimental data-.on the triangular lattice that are most prevalent at 1/2, followed by 1/4, 1/3, 3/8, 
and so on. The experimental datatl also shows that the system has four different regions of superconducting stability 
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for one value of t, centered at 0, 1/2, 1/4 and 1/3. While our results in Fig. ||(a) only show superconductivity around 
4> = [where A(t) = 0], we can see that if t was increased, and the calculations carried out to higher order, it is 
likely that we would see additional superconducting regions appearing (we believe first around 1/2, followed by 1/4 
and 1/3). 

Finally we study the magnetic-field dependence of the dynamical critical exponent zv (Eq. [?]) for small fields in 
Fig. ||(b). The BHM in zero magnetic field can be mapped onto a three-dimensional XY model which has zv = 0.67 
independent of the lattice structurelI3. By using the Pade analysis on our third-order expansion, we obtain zv = 0.6-L 
for the square lattice and zv = 0.72 for the triangular lattice at 4> = 0. A Pade analysis of a tenth-order expansions 
yields zv = 0.69 for both lattice types, which shows the convergence of higher-order calculations in our method. As 
seen in Fig. |[ the dynamical exponent appears to increase as the magnetic field increases. It is, however, difficult to 
conclude whether zv remains equal to it's zero-field value, increases with magnetic field, or immediately jumps to one 
upon the onset of a magnetic field, solely on the basis of this third-order analysis. But, the fact that the most likely 
point to have zv < 1 is <j) = 1/2 (because it requires only two components for the order parameter), and a higher-order 
expansion predicts zv — 1 there, leads us to conjecture that zv — 1 for all nonzero magnetic fields. This latter result 
is supported by Monte Carlo data on the antiferromafflaetic XY model for stacked triangular planes, which have a 
weakly first-order transition rather than a critical pohralS. 

IV. CONCLUSIONS 

In conclusion, we applied a strong-coupling i/L7-expansion (up to third order) to study the insulator-superfluid 
phase transitions of the two-dimensional Bosc-Hubbard model under the influence of a magnetic field. Although 
our analysis is limited, we are able to produce reasonable results that both agree with physical intuition and with 
experiments. We found that the Mott insulating phase enlarges with an increasing magnetic field. This is explained 
by the localizing effect of the magnetic field on the moving bosons. Qualitative agreement in the increase of the critical 
hopping energy t C rit was found with experimental results on Josephson junction arraysa. For small magnetic fields, we 
find a power-law behavior of the gap energy close to the critical point (t crit , pL crit ) with a dynamical critical exponent 
that increases with <f>. For larger magnetic fields, we find a repulsion of the two phase boundaries, which indicates 
a change from the power-law like behavior to either Kosterlitz-Thouless behavior or to a "first-order" transition. 
Our results are also consistent with the "critical point" immediately changing to a "first-order" discontinuous change 
in the slope of the phase boundaries as the magnetic field is turned on. We found that the gap energies for small 
fixed hopping and variable magnetic field illustrate commensurability effects for rational fluxes which is also seen in 
Josephson-junction array experiments. More work is needed to understand the change in character of the insulator 
to superfluid phase transition as a magnetic field is introduced: does the system have a field-dependent power-law 
dependence which crosses over to a "first-order" (or more exotic cusp-like shape) as the field increases, or does it 
immediately become "first-order" in a field. Higher-order calculations are needed to decide this issue. It would also 
be interesting to extend the scaling analysis of the BHM to include its behavior in an external magnetic field. 
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